Mon. Not. R. Astron. Soc. 001, 1-17 (2012) Printed 24 February 2013 (MN KTeX style file v2.2) 

Prospects of observing a quasar HII region during the Epoch of 
Reionization with redshifted 21cm 



Kanan K. Datta^ *, Martina M. Friedrich\ Garrelt Mellema\ Ilian T. Iliev^, Paul R. Shapiro^ 

^ Department of Astronomy & Oskar Klein Centre, AlbaNova, Stockholm University, SE-106 91 Stockholm, Sweden 

^Astronomy Centre, Department of Physics & Astronomy, Pevensey II Building, University of Sussex, Falmer, Brighton BNl 9QH 

^Department of Astronomy and Texas Cosmology Center, University of Texas, Austin, TX 78712, USA 



(N 

O 
<N 



24 February 2013 



o 
u 

o 

<N 
> 

CO 

O 
(N 



^ 



ABSTRACT 

We present a study of the impact of a bright quasar on the redshifted 21cm signal during the 
Epoch of Reionization (EoR). Using three different cosmological radiative transfer simula- 
tions, we investigate if quasars are capable of substantially changing the size and morphology 
of the H II regions they are bom in. We choose stellar and quasar luminosities in a way that is 
favourable to seeing such an effect. We find that even the most luminous of our quasar models 
is not able to increase the size of its native H II region substantially beyond those of large 
H II regions produced by clustered stellar sources alone. However, the quasar H II region is 
found to be more spherical. We next investigate the prospects of detecting such H II regions 
in the redshifted 21cm data from the Low Frequency Array (LOFAR) by means of a matched 
filter technique. We find that H II regions with radii ~ 25 comoving Mpc or larger should 
have a sufficiently high detection probability for 1200 hours of integration time. Although the 
matched filter can in principle distinguish between more and less spherical regions, we find 
that when including realistic system noise this distinction can no longer be made. The strong 
foregrounds are found not to pose a problem for the matched filter technique. We also demon- 
strate that when the quasar position is known, the redshifted 21cm data can still be used to 
set upper limits on the ionizing photon rate of the quasar. If both the quasar position and its 
luminosity are known, the redshifted 21 cm data can set new constrains on quasar lifetimes. 
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1 INTRODUCTION 

The Epoch of Reionization (EoR) refers to the timespan between: 
(1) the formation of the first sources of light in the first halos that 
collapse under the influence of self-gravity in the highest density 
peaks in the neutral Universe after recombination; and (2) the time 
when most of the intergalactic hydrogen has become ionized by the 
radiation emitted from these sources of light. There are currently 
two key observations constraining this timespan. The optical depth, 
Tes, due to Thompson scattering off free electrons, is measured by 
means of the polarization power spectrum of the cosmic microwave 
background (CMB). These measurements have constrained res to 
0.088 ±0.015, implying that an instantaneous reionization would 
have happened at ^ = 10.4 ±1.2 (Komatsu et al. 2010). Quasar 
spectra obtained within the Sloan Digital Sky Survey (SDSS) indi- 
cate a low, but rapidly rising neutral fraction around redshift 6 (Fan 
et al. 2006; Willott et al. 2007). The combination of these two mea- 
surements suggests that the epoch of reionization extended over at 
least several redshift units. 
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Radio telescopes capable of measuring at low frequencies, 
(GMRT\ 21CMA^ LOFAR^ MWA^, PAPER^) are either already 
existent or in the process of being build and aim to detect the signa- 
ture of redshifted 21cm radiation from neutral hydrogen during the 
EoR. These measurements will provide completely new constraints 
on this early epoch of galaxy formation. 

The sources of the ionizing photons during the EoR are po- 
tentially a mixture of first galaxies, metal free POP III stars and 
quasars. Although the contribution from quasars is believed to be 
small (see e.g. Giroux & Shapiro 1996; Loeb 2009) due to their de- 
creasing space density above redshift 2-3 (e.g. Schmidt et al. 1995; 
Boyle et al. 2000; Cristiani et al. 2004; Gonzalez- Serrano et al. 
2005), there are still doubts about the role of quasars in the EoR 
since secondary ionizations can yield several tens of ionizations 
per emitted high energy photon (Volonteri & Gnedin 2009) and 



Giant Metrewave Telescope, http://gmrt.ncra.tifr.res.in 

21 Centimeter Array, http://21cma.bao.ac.cn 

Low Frequency Array, http://www.lofar.org 

Murchison Widefield Array, http://www.mwatelescope.org 

Precision Array to Probe the EoR, http://astro.berkeley.edu/~dbacker/eor 



© 2012 RAS 



Datta et al. 



because galaxies alone might not be sufficient suppliers of ioniz- 
ing photons (Bunker et al. 2010). Srbinovsky & Wyithe (2007) use 
a model for the quasar luminosity function together with a semi- 
analytic reionization description together with observations to set 
an upper limit to the quasar contribution to the ionizing background 
at 2; ^ 5.8 and find an upper limit of 14%. Based on quasar survey 
data, Willott et al. (2010) estimated the photon rate density from 
quasars to be between 20 and 100 times lower than the required 
rate to complete the reionization. Dijkstra et al. (2004) and Dijkstra 

(2006) claim that the observed X-ray background can be used to 
set a limit on the contribution of power law type sources to global 
reionization. 

Apart from the observational limits on the number of quasars, 
it also remains unclear whether sufficiently massive black holes can 
form rapidly enough in order to create luminous quasars at redshifts 
relevant for hydrogen reionization. Volonteri (2010) suggests sev- 
eral different formation mechanisms for massive black holes. The 
recent discovery of a very luminous quasar at 2; ^-^ 7 (Mortlock et al. 
201 1) shows however that there are in fact very massive black holes 
at such high redshifts. Mortlock et al. (2011) estimate the mass of 
the central black hole to be of the order of 10^ M© based on Mg II 
line width measurements and the luminosity at A = 300 nm using 
the method described in Vestergaard & Osmer (2009). 

This highest redshift quasar was detected as part of the UKIRT 
Infrared Deep Sky Survey (UKIDSS) which is expected to find 
'-^ 10 more quasars dX z > 6.5 (Lawrence et al. 2007). Further 
surveys for high redshift quasars are being performed with ESO's 
VISTA telescope (Emerson & Sutherland 2010), the Panoramic 
Survey Telescope and Rapid Response System 1 (PAN-STARRS 
Kaiser et al. 2002; Morganson et al. 2011) and the Canadian- 
France-Hawaii Telescope (the CFHQSIR program). These surveys 
will add to the already existing list of high redshift quasars estab- 
Ushed by the Sloan Digital Sky Survey (SDSS Fan et al. 2001, 
2006; Jiang et al. 2008) and the Canadian-France High Redshift 
QSO survey (CFHQS Willott et al. 2007, 2010). All of the above 
efforts concentrate on the redshift interval 6 to 7. To map out the 
quasar population for 2; < 8 a space mission such as ESA's EU- 
CLID may be needed (Roche et al. 2012). 

While the quasar contribution to the ionizing photon budget 
is a matter of ongoing debate, it is generally accepted that quasars 
(including the so called mini- and micro-quasars) contribute sub- 
stantially to the (global) heating of the neutral IGM before large 
scale ionization. Several EoR simulations (e.g. Back et al. 2010; 
Zaroubi et al. 2007) include a quasar-like population to investigate 
effects of heating. There are other groups that include a mixture of 
sources to test the large scale effect of quasars on the morphology 
of H II regions, for example Geil & Wyithe (2009) who include a 
quasar population in a semi-analytic fashion and investigate their 
effect on the 21cm power spectrum. 

A different approach is taken by Maselli et al. (2007), who 
include a single quasar in a reionization simulation to generate 
Lya absorption spectra. They concentrate on the very late stages of 
reionization, around z > 6, and assume an IGM in photoionization 
equilibrium with some uniform ionizing background. Their main 
goal is to investigate the reliability of measurements of the average 
ionization state of the IGM from estimates of near-zone sizes from 
quasar spectra. They come to a similar conclusion as Lidz et al. 

(2007) who use a semi-numerical approach to model quasar spec- 
tra from a quasar in a patchy ionized IGM, namely that the effect 
of small scale inhomogeneities in the IGM is larger than the effect 
of the global ionization fraction and therefore conclude that Lya 
absorption spectra of high redshift quasars give only limited infor- 



mation about the global ionization fraction. Additionally, Alvarez 
& Abel (2007) point out the importance of the pre-existence of stel- 
lar H II regions at the locations of quasars when they begin to shine. 
Neglecting the stellar H II region will will lead to an overestimate 
of the calculated average surrounding ionization fraction. 

A similar approach of concentrating on a single quasar has 
been taken recently by Majumdar et al. (201 1) and earlier for exam- 
ple by Geil & Wyithe (2008). While the latter investigate the effect 
of quasar H II regions on mock 21cm line-of-sight spectra to extract 
information about the global ionization fraction, the former use a 
filter technique described first in Datta et al. (2007) to investigate 
quasar properties from 21cm observations. Our paper has a similar 
goal, namely to study the imprint of a single bright quasar on the 
morphology of H II regions and its observability using 21cm ob- 
servations. The main difference between these previous approaches 
and ours is that they use a semi-numerical method for generating 
the ionization fraction field including the quasar while we perform 
full radiative transfer in a hydrogen/helium IGM. 

This paper tries to answer the following questions: Does a 
quasar change the size and the morphology of an H II region in 
a distinct manner? Can LOFAR in combination with the matched 
filter approach be used to detect individual H II regions? Can we 
use the matched filter technique in a blind search (i.e. without prior 
knowledge of a quasar location) to identify an H II region as a 
quasar H II region? If we have prior knowledge about the quasar 
position, can we use the matched filter technique to extract quasar 
properties from the H II region size estimated? To answer these 
questions, we test three cases of detectability of a single quasar 
H II region in a patchy ionized environment. 

The structure of the paper is as follows. In Sect. 2 we de- 
scribe our simulations. In Sect. 3 we compare the ionization frac- 
tion fields and differential brightness temperature maps of the sim- 
ulations without a quasar to the ones with a quasar and make some 
estimate on the expected quasar H II region sizes based on photon 
counts. We try to asses the three dimensional shape of the H II re- 
gions. Sect. 4 describes how to simulate the radio-interferometric 
visibilities from the redshifted 21cm line radiation. Sect. 5 summa- 
rizes the matched filter method which we use to detect the quasar 
H II regions in Sect. 6. Sect. 7 discusses the role of foregrounds. In 
Sect. 8 we summarize our findings. 



2 SIMULATIONS 

The general outline of our simulations that we perform to test the 
detectability of quasar H II regions is the following: 

(a) We start with the results of an already existing large scale hydro- 
gen only radiative transfer reionization simulation with stellar-type 
sources. 

(b) In this simulation, we select a redshift Zon where the global ion- 
ization fraction is small enough that the individual H II regions of 
galaxy clusters, although already partly connecting, can still be rec- 
ognized as individual H II regions. 

(c) We add helium, initialized to have the same ionization state as 
hydrogen some millions of years before the quasar turn on and run 
the simulation until we reach Zon in order to have better initial con- 
ditions for helium when including the quasar. 

(d) At Zon, we select the most massive halo. Here we include one 
additional source, a quasar with an ionizing photon rate that is con- 
stant in time and which is chosen on the basis of the host-halo mass. 

(e) We run the simulation for another 23 Myr with all stellar sources 



© 2012 RAS, MNRAS 001, 1-17 



Prospects of observing a quasar HII region during the EoR region with redshifted 21cm 3 



and the quasar, following the ionization of hydrogen and helium. 
In the following, we describe these steps in more detail. 



2.1 Simulations without a quasar 

The simulation methodology (with stellar type sources only) has 
been described in detail in Iliev et al. (2006), Mellema et al. (2006a) 
and Iliev et al. (2007). Here we only summarize some important as- 
pects. The density fields through which we trace the ionizing radi- 
ation and the sources that emit this ionizing radiation are extracted 
from cosmological N-body simulations of structure formation. For 
these N-body simulation, the CubeP^M code was used. This code 
was developed from the PMFAST code (Merz et al. 2005), see Iliev 
et al. (2008b) for a short description of the CubeP^M code. It uses 
particle-particle interactions at sub grid distances and particle-mesh 
interactions for larger distances. We use the results of two such 
simulations: one with a volume of (163 cMpc)^ that has 3072^ 
particles and a mesh size of 6144^ cells; and one with a vol- 
ume of (607 cMpc)^ that has 5488^ particles and a mesh size of 
10976^ cells.^ This implies particle masses of 5.5 x 10^ M© and 
5.0 X 10^ M© , respectively. The first set of parameters guarantees a 
minimum resolved halo mass of ^^ 10^ M© which is approximately 
the minimum mass of halos able to cool by atomic hydrogen cool- 
ing. The second set of parameters guarantees a minimum resolved 
halo mass of ^ 10^ M©. 

All N-body simulations use the cosmological parame- 
ters for a flat ACDM universe with (Qrn, ^b, /i,^, crs) = 
(0.27, 0.044, 0.7, 0.96, 0.8), based on the five year WMAP results 
(Komatsu et al. 2009). For creating the halo list, a spherical over- 
density halo finder was used. 

The density is assigned to a coarser mesh (256^ and 504^, re- 
spectively) by smoothing the dark matter particle distribution using 
an SPH kernel function. To convert the dark matter density into 
a baryon density, we assume that the gas distribution follows the 
dark matter. This is valid on the scales of the radiative transfer cells 
(0.64 and 1.2 comoving Mpc, respectively) since the local Jeans 
length at the mean density of the IGM is much smaller than that. 
For the radiative transfer simulations, we use the short character- 
istic ray tracing method C^-Ray (Mellema et al. 2006a; Friedrich 
et al. 2012) on this coarser mesh. 

Each halo of mass M is assigned an ionizing photon luminos- 
ity 

MQb 



N^ =^7 



lOMyrQmnip ' 



(1) 



where N^ is the number of ionizing photons emitted per Myr, M 
is the halo mass, and rrip is the proton mass. Halos are assigned 
different luminosities according to whether their mass is above or 
below 10^ M© (but above 10^ M©). For the second simulation, we 
use a recipe calibrated against the 163 cMpc simulation, to include 
sources in halos with masses lO^M© < M < 10^ M©. Sources in 
low mass halos are suppressed if the ionization fraction of their cell 
is larger than 10%. In the simulations we use here, sources in high 
mass halos have an efficiency g^^ — 1.7 and sources in low mass 
halos have an efficiency g^ — 8.7. 

These are rather low efficiencies for the stellar sources. This 
combined with our prescription for the quasar luminosity (see Sect. 
2.2) means that the quasar will outshine the stellar sources in its 



^ As is described in detail below, we only resimulate a sub-volume of the 
latter. 



Table 1. Summary of important quantities of the three simulated cases. 

resolution box (sub-box) ^(QSO on) ^(QSO on 2;(QS0 on 

size +1 1 .5 Myr) +23.0 Myr) 

early QSO 0.64 cMpc 163 cMpc 8.636 8.515 8.397 

lateQSO 0.64 cMpc 163 cMpc 7.760 7.664 7.570 

large box 1.2cMpc 607 (242) cMpc 7.760 7.664 7.570 




Shanks et al. 201 1 


O X=1 , early QSO turn on 


X=0.333, early QSO turn on 


X ?i=1, late QSO turn on 
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— Ferrarese (2002), Eq. 7 


Ferrarese (2002), Eq. 6 


Ferrarese (2002), Eq. 4 



Figure 1. Different analytic estimates for the halo mass M^ - black hole 
mass Mbh relation. Note the large spread. The Shankar et al. (2010) 
curves are produced using their fiducial parameters. Shanks et al. (2011) 
and Croom et al. (2005) suggest a narrow host halo mass range for quasars, 
indicated in the figure by the grey rectangular shaped region. Included in 
this plot are the quasars from our simulations. We include two possible 
black hole masses for each quasar implementation, one assuming A = 1 
and one assuming A = 0.333, which both produce the same luminosity. 
DeGraf et al. (2012) find that the distribution of Eddington ratios peaks 
between A = 0.1 — 0.2 for black hole masses between Mbh = (10^- 
— 1O^)M0. The halo mass for the six points corresponds to the mass of 
the most massive halo for the indicated cases. 



host halo by a factor of a few, see Table 2 for a summary of stellar 
and quasar photon rates. This is necessary for the quasar to have 
any impact on the neutral hydrogen distribution in its environment 
as otherwise the stellar sources in neighbouring halos would always 
dominate the photon budget in the H II region surrounding this high 
density region. 

The first simulation used here, 163Mpc_gl.7_8.7S is one sim- 
ulation of a set of simulations presented and analysed in Iliev et al. 
(2011). It has a simulation volume of (163 cMpc)^. The second 
simulation, 607Mpc_gl.7_8.7S has a considerably larger volume. It 
will be analysed and described in detail in Iliev et al., in prep. 

The 163 cMpc hydrogen only simulation without the QSO 
reaches global ionization fractions (xm)=(0.1, 0.5, 0.7, 0.9 and 
0.99) at redshifts z =(9.9, 7.6, 7.2, 6.9 and 6.7) (see Iliev et al. 
2011). The 607 cMpc simulation has currently not reached more 
than (xm) =0.7, the redshifts for global ionization fractions (0.1, 
0.5, 0.7) are (9.6, 7.5, 7.1). Although the two simulations were per- 
formed using the same source properties, it is visible that the ion- 
ization history is slightly delayed in the larger simulation volume. 
This is most probably due to the way in which the suppressible 
sources are included, which results globally in a higher suppressed 
fraction. The 163 cMpc simulation has an electron Thompson scat- 
tering optical depth r = 0.058. 



2.2 Adding a quasar to the simulations 

For adding a quasar, the next step to take is to identify a suitable 
redshift. Observationally, the highest redshift quasar has been de- 
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tected dii z — 7.085 (Mortlock et al. 2011). However, we choose 
higher redshifts {z = 8.636 and z — 7.760 for the quasar turn on) 
because of the lower global (i.e. mass averaged) ionized fractions 
(xm) at those redshifts (0.25 and 0.45) in our simulation. At these 
epochs the size distribution of H II regions still allows them to be 
rather isolated and thus makes it possible for a quasar to modify 
the H II region that it finds itself in. As was shown for example 
by Friedrich et al. (2011) at higher global ionized fractions most 
H II regions have become connected and the topology of the IGM 
is rather complex. Geil & Wyithe (2008) also already concluded 
that beyond (xm) ^ 0.7 the signature of the quasar is difficult to 
recognize in the neutral hydrogen distribution. 

As was explained in Friedrich et al. (201 1) a reionization sim- 
ulation of a given volume can only accommodate a limited number 
of H II regions of a given size at a given global ionized fraction. 
This means for our case that the (163 cMpc)^ volume may possi- 
bly only accommodate one large H II region from which we would 
then conclude that the H II region containing the quasar could be 
recognized as being the largest one. We therefore also use the larger 
volume of (607 cMpc)^ which allows us to address the question 
whether there are H II regions powered by stellar sources alone 
which could obtain a similar size as the quasar containing H II re- 
gion. An additional advantage of using this larger volume is that 
the most massive halo corresponds to a higher cr-peak in the den- 
sity field and is therefore more massive. 

However, because of the larger volume, our procedure for 
adding the quasar to the 607 cMpc box differs a little from that for 
the 163 cMpc one. Firstly, the radiative transfer of ionizing photons 
in the larger volume is limited by a mean free path set to 60 cMpc, 
roughly corresponding to the mean free path of ionizing photons 
due to Lyman Limit Systems dX z ^ Q (Songaila & Cowie 2010). 
We therefore only need to resimulate part of the volume to assess 
the impact of the quasar as photons further away than 60 cMpc 
from the region of interest will not reach this region. In order to in- 
sert the quasar region back into the whole volume we actually need 
to include a buffer zone of the same extend around our resimulated 
region, so the "zoomed" region which we resimulate with a quasar 
has a size of (242 cMpc)^ centred on the quasar and consists of 
201^ radiative transfer cells. This procedure introduces the follow- 
ing differences compared to the quasar simulations in the smaller 
volume: 

(1) By taking a sub-box, the average density in the box is no longer 
equal to the average density of the universe. In fact, with our choice 
of size, the sub-box is slightly less dense than the average density 
of the universe. However, this is only a minor deviation, the average 
density in the sub-box is only 0.25% lower than the average density 
in the whole simulation volume. 

(2) The spatial resolution of this simulation is worse than for the 
(163 cMpc)^ simulation, 1.2 cMpc instead of 0.64 cMpc. 

(3) Since we take a sub- volume of a larger simulation, we use open 
boundary conditions. When inserting the sub volume back into the 
full volume we discard the outer buffer zone. 

The basic simulations described above, 163Mpc_gl.7_8.7S 
and 607Mpc_gl.7_8.7S do not include helium. This is equivalent 
to assuming that helium is singly ionized everywhere where hy- 
drogen is ionized and that double ionization of helium is neg- 
ligible. This is an assumption frequently made in EoR simula- 
tions including stellar type sources (e.g. Iliev et al. 2006; Mellema 
et al. 2006b; McQuinn et al. 2007; Zahn et al. 2007; Aubert & 
Teyssier 2010). As was demonstrated in Friedrich et al. (2012), 
this is a valid assumption in the case of stellar type sources. Since 
quasars emit a substantial amount of their photon output at en- 



ergies that are sufficient to doubly ionize helium, and since the 
ionization cross-section of helium is roughly an order of magni- 
tude larger than that of hydrogen, it is essential to include helium 
in our simulations containing a quasar. In order to do so, we ini- 
tialize the helium ionization fractions 11.5 Myr before the quasar 
turn on for each of the three cases to correspond to that of hy- 
drogen. He III is assumed to be absent. We showed in Friedrich 
et al. (2012) that the He I fraction field corresponds approxi- 
mately to the H I fraction field. We evolve the so generated he- 
lium ionization fraction field together with the hydrogen for 11.5 
Myr before adding the quasar source to get better initial condi- 
tions for the helium fraction fields. The corresponding redshifts 
for adding helium are 8.762 for the early quasar turn of simulation 
163Mpc_gl.7_8.7S and 7.859 for the late quasar turn on of simula- 
tion 163Mpc_gl.7_8.7S and simulation 607Mpc_gl.7_8.7S. In Ta- 
ble 1 we summarize all relevant quantities for all three (early and 
late QSO turn on of simulation 163Mpc_gl.7_8.7S and the large 
box simulation 607Mpc_gl.7_8.7S) studied cases. 

The general idea for the quasar model is to connect the quasar 
luminosity to the host halo mass. This can be done in three steps: 

(1) The quasar luminosity Lqso can be connected to the black hole 
mass Mbh (reverberation mapping shows that there is a correlation 
between black hole mass and quasar luminosity, Shen et al. 2008) 

(2) Mbh can be connected to the mass of galaxy bulges Mb. (Mag- 
gorian type relation, Magorrian et al. 1998; Marconi & Hunt 2003; 
Haring & Rix 2004; Tundo et al. 2007; Graham 2012). 

(3) Mb can be connected to the total mass of the halo Mh (Fer- 
rarese 2002). 

Combining these relations, there is a big scatter for the relation 
halo mass to black hole mass, see Fig. 1. Shanks et al. (2011) ar- 
gue that such a general relation does not exist. Instead, they suggest 
that there is a typical halo mass range for active quasars. Also Kor- 
mendy & Bender (201 1) argue that there is in general no correlation 
between the halo mass of arbitrary galaxies and the central black 
hole. However, due to the lack of better observational or theoretical 
restrictions, we base our quasar luminosity on the host halo mass. 

Since we aim for a luminous quasar, we orientate us towards 
the upper limit of the conversion factor between halo mass and 
black hole mass, see Fig. 1. The mass of the black hole can be 
converted into a luminosity by assuming an Eddington limited ac- 
creation 

Lqso = AL,:,, = a '"^^""^^^ « a 1.3 X 10^^ ^W. (2) 



(7T 



Mrr 



Here, G is the gravitational constant, Mbh the mass of the black 
hole, TTip the proton mass, c the speed of light and ax the Thomp- 
son scattering cross section. A is the Eddington efficiency parame- 
ter. Although Eddington efficiencies A larger than 1 are in principle 
possible (see for example Heinzeller et al. 2007, for simulations as- 
suming a thin disc), observations show that most quasars are below 
the Eddington limit (see for example figure 1 in Steinhardt & Elvis 
2010). Also hydrodynamical simulations indicate that most quasars 
shine with luminosities below the Eddington luminosity (DeGraf 
et al. 2012). In Fig. 1, we include the corresponding black hole 
masses for two assumed values for A, to give an impression of how 
our choice of parameters compares with other models. The impor- 
tant parameter for us is the total luminosity and not the black hole 
mass, so there is a degeneracy between A and Mbh- Therefore, 
the six points in the plot only correspond to three different sim- 
ulations with three different quasar luminosities. The conversion 
factor between halo mass and quasar luminosity, Lqso = ACMh 
for the early quasar turn on is AC = 4.3 x IO^^WMq^ and 



© 2012 RAS, MNRAS 001, 1-17 



Prospects of observing a quasar HII region during the EoR region with redshifted 21cm 



Table 2. Summary of important results for the three cases; The stellar ion- 
izing photon rates { A^* } are averages over the stellar photons emitted in the 
quasar area of influence during the quasar-on time, A^* and N^^^ denote 
the total number of stellar and quasar photons emitted in the region of in- 
fluence of the quasar over its entire history, respectively. M^ ^lax denotes 
the mass of the most massive halo in each simulation volume. 



(iV;)/s 



iV:^'°/s Mh,max/M0 



NX 



NT° 



early QSO 4.5 x 10^^ 3.3x10^5 4.9x10^1 2.5 x 10"^° 2.5x10^° 
lateQSO 1.7 x lO^^ 1.4x10^6 6.9x10^1 9.8 x 10"^° 1.0 x lO'^^ 
large box 4.3x10^5 2.4x10^6 1.2x10^2 2.2 x lO"^! 1.7 x 10^^ 



for the late quasar turn on in the 163 cMpc and 607 cMpc sim- 



1.3 X 10^^WM~. The most massive halo in 



ulations AC 

the three cases (early turn on, late turn on and large box) are 
4.9 X lO^^Mo, 6.9 X lO^^M© and 1.2 x lO^^M©, respectively. 
With our choice of power law index (a = — 1.5, L(i/) oc v^) 
and with our high energy cut off at 5440 eV, this corresponds to 
3.3 X 10^^ s"\ 1.4 X 10^^ s"^ and 2.4 x 10^^ s"^ ionizing 
photons per second for the early and late quasar turn on of the 
163 cMpc simulations and the 607 cMpc simulation, respectively. 
We summarize these numbers in Table 2. The recently observed 
z = 7.085 quasar is estimated to have an ionizing photon rate of 
1.3 X 10^^ s-^ (Mortlock et al. 2011). The lower redshift (z - 6) 
quasar that Maselli et al. (2007) investigate in their simulations has 
an ionizing photon rate of 5.2 x 10^^ s~^ and the quasar included 
in the simulations of Majumdar et al. (201 1) has an ionizing photon 
rate of 8.97 x 10^^ s~^. Although our most luminous quasar has a 
lower ionizing photon rate by a factor of a few compared to these 
other studies, the resulting H II regions have comparable sizes since 
we consider a longer quasar lifetime. The reason we have a slightly 
less luminous quasar is that our most massive halo in our simu- 
lation is less massive than the most massive halo in, for example, 
Maselli et al. (2007) since we are considering higher redshifts. 

Cano-Diaz et al. (2012) observe a 2; = 2.4 quasar and find that 
star formation in the host halo is quenched. We do not include neg- 
ative feedback from the quasar on the star formation other than the 
usual suppression of low-mass sources in ionized regions. How- 
ever, since the quasar dominates the photon output during its on- 
time, see Table 2, and since Iliev et al. (2008a) found that the con- 
tribution of the largest halo to the total photon-budget in a jointly 
produced H II region is small compared to the sum of photons from 
the clustered fainter sources contributing to the H II region, we ar- 
gue that negative feedback on the star formation would not change 
our conclusions significantly. Another important assumption is that 
the quasar emits radiation equally in all directions. Although this is 
in line with previous studies of this kind (Geil & Wyithe 2008; Ma- 
jumdar et al. 201 1), real quasars may emit non-isotropically, which 
would change the geometry of their region of influence. Consider- 
ing this complication is beyond the scope of this work. 



3 RESULTS FROM THE SIMULATIONS 

In this section, we first investigate the effect of the quasar on the 
size and morphology of the H II region it is born in. We estimate 
the expected size of the quasar H II region on the basis of photon 
counts. We compare the morphology of the H II region around the 
most massive halo with and without the quasar included by means 
of the ionization fraction fields and column density maps. Next, we 
compare the quasar H II region for one of the three cases to a large 
H II region solely made by stellar sources. We show 21cm maps 
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Figure 3. Column densities along the line of sight over a distance roughly 
corresponding to the expected joint H II region diameter as indicated in the 
plot. Panel (a) shows the early quasar turn on case, (b) the late quasar turn 
on case and (c) the large box case. The green squares in panel (c) indicate 
the size of the smaller simulation volume (163 cMpc). 



with a smoothing that corresponds roughly to LOFAR resolution. 
To quantify the changes in size and morphology, we use a statis- 
tical measure using many lines of sight from the central source, 
introduced by Iliev et al. (2008a). Finally we analyse the He III 
region around the quasar. 
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Figure 2. Slices of hydrogen ionization fraction, blue corresponds to fully neutral, red to fully ionized. The upper left four slices (panel a) correspond to the 
early quasar turn on in the 163Mpc_gl.7_8.7S simulation, the upper right four slices (panel b) correspond to the late quasar turn on in this simulation. The 
lower right four slices (panel d) correspond to a sub-region of the same size (163 Mpc) of the 607Mpc_gl.7_8.7S simulation. In each panel, the upper row 
correspond to a simulation without quasar and the lower row to the same simulation including a quasar in the centre. For completeness, we also show a slice 
of the whole 607Mpc_gl.7_8.7S (panel c) simulation. The yellow square indicates the re-simulated part and the green square indicates the sub-region show in 
panel d. 



3.1 



Expected sizes of quasar H II regions from photon 



The emitted photon rate of the quasar for each case, see Table 2, 
can be used to estimate very roughly the comoving radius R of 
an expected spherical H II region after time ton if the surrounding 
medium is assumed to be completely neutral: 



R 



AtttIc 



counts where ton is the time since the quasar turned on and ric is the av- 

erage comoving number density (hydrogen+helium) of the IGM. 
Evaluating this for the early and late quasar turn on of simulation 
163Mpc_gl.7_8.7S and for the large box using ton = 23.0 Myr, 
gives roughly 10 cMpc, 16 cMpc and 19 cMpc, respectively. 

However, all stellar sources in this region will also contribute 
to the "quasar" H II region. For the early quasar turn on, the num- 
ber of stellar photons in this region from the formation of the first 
source in this region to the end of the quasar phase is: N* ^ 
2.5 X 10^*^. The total number of emitted ionizing photons from the 
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quasar is N^^° ^ 2A x 10^°. In analogy, for the late quasar turn 
on the numbers are A^* - 9.8 x 10'^° and Ar^"° - 1.0 x 10'^\ The 
large box simulation has N^ - 2.2 x 10^^ and A^^^ - 1.7 x 10^\ 
These numbers are summarized in Table 2 

We note that the contribution from the quasar in terms of total 
number of emitted photons in the region of influence of the quasar 
is about 50%. We also note that the ratio A^^ /N^^"" is largest for the 
most massive of our quasars. As Geil & Wyithe (2008) mention, 
the effect of suppression of low mass sources delaying reionization 
globally, can be reduced in highly overdense regions due to the bi- 
ased formation of very massive galaxies which are not suppressed. 

Using the total number of photons, i.e. N^ + N^^^ , for es- 
timating the comoving radius for the joined H II region gives 
R ^ 12.0 cMpc for the early quasar turn on and R ^ 19.7 cMpc 
for the late quasar turn on and 25.1 cMpc for the large box. ^ 



3.2 Comparisons between the morphology of H II regions 
with and without quasar 

In Fig.2, we show for each of the three cases (early quasar turn on 
in panel (a), late quasar turn on in panel (b), large box in panel (d)) a 
slice in the plane of the most massive halo hosting the quasar at 1 1 .5 
Myr and 23.0 Myr after the quasar turn on (left and right columns 
of each panel, respectively). For comparison, for each case at the 
same times, we also show that slice in the simulation without the 
quasar (upper row in each panel). For better comparison, we ad- 
justed the size of the slices from the 607Mpc_gl.7_8.7S simulation 
to a side length of 163 cMpc centred around the quasar. In panel 
(c), we show the full size of the (607 cMpc)^ box and indicate the 
re- simulated region and the sub-region shown in panel (d). 

It is important to note that we show only one 2-dimensional 
slice through the ionization fraction field for each case. From this 
representation alone, no absolute judgement about the size, spheric- 
ity or connectivity of the H II region should be made. However, we 
used among other things iso- surface representations and the above 
mentioned line-of-sight statistic (described in detail below) to con- 
firm the statements below. It is useful to note the following from 
the simulations without a quasar (upper rows of panel a, b, d and 
panel c): The H II region around the most massive halo is large 
compared to average H II regions but not necessarily the largest 
one. Especially in the large box (panel c), there is an H II region 
(marked with "B") which is larger. The reason for this is that there 
are several halos in the simulation volume which are competing for 
the most massive halo between 2; ^-^ 30 — 6. They are located in 
equally overdense regions with many clustered sources. 

From the ionization fraction fields with quasar (lower rows in 
panels (a), (b) and (d)), we note that the H II regions of the quasars 
have as sharp boundaries as the stellar H II regions. This is con- 
sistent with previous results by Shapiro et al. (2004), Thomas & 
Zaroubi (2008) and the more detailed study of Kramer & Haiman 
(2008) who found that a harder spectrum (a > —1.2) and some 
level of obscuration of the low frequency quasar photons are needed 
to result in a detectable ionization front thickness. 

Furthermore, it is visible that the growth of the H II region 
during the quasar lifetime is dominated by the quasar (by construc- 
tion). As we showed above, the total photon budget from the time 
of the formation of the first source in that region until the end of the 



^ Since this corresponds to approximately 3, 8 and 9.5 x 10^ proper lyr, 
respectively, the finite speed of light is not an issue since this is below the 
timescale we are interested in. 



quasar lifetime is not dominated by the quasar. We summarize the 
numbers is Table 2. 

Next, we consider the neutral hydrogen density. We calculate 
the neutral column density along the line of sight for a distance 
corresponding to the expected H II region diameter (from photon 
counts, see above) centred around the most massive halo, see Fig. 
3. Assuming the quasar to be observed at 2; = 8.397 and z — 7.570 
respectively, the distance over which we calculated the H I column 
density corresponds to Az = 0.082, 0.125 and 0.150 for the early, 
late and large case, respectively. This corresponds to the following 
frequency ranges. For the early case (150.451 - 151.769) MHz; for 
the late case (164.532 - 166.949) MHz and for the large box case 
(164.285 - 167.185) MHz. This corresponds to an integration over 
roughly 1300, 2400 and 2900 intrinsic LOFAR frequency channels 
of 1 kHz. However, due to storage restrictions, LOFAR analysis 
will probably use channels of 100 - 500 kHz. 

The size along the line of sight also implies a time difference 
from the far side to the quasar of about 3 Myr, 8 Myr and 9.5 Myr, 
for the early, late and large box case, respectively. Therefore, we ac- 
tually see the far-side of the quasar H II region at a stage that corre- 
sponds more to the one shown in the left columns of panels (a), (b) 
and (d) in Fig. 2. However, we neglect this effect and consider the 
quasar H II region everywhere after 23.0 Myr. The main reason for 
this is that our time resolution is not sufficient to make meaningful 
interpolations. Assuming the quasar turns off at the time of obser- 
vation, the side of the H II regions closer to us starts recombining. 
However, the recombination time, even in overdense regions with 
5 r^ 10 is still of the order of 50 Myr. For comparison, the overden- 
sity in the H II region is less than 0. 1 in the large box case. We note 
that with our particular choice of lifetime and photon output, the 
part of the quasar H II region furthest from us had already time to 
grow to a non-negligible size (see left columns, lower row of pan- 
els (a), (b) and (d) in Fig. 2). Also, less dependent on the particular 
quasar lifetime, the pre-existing stellar H II region itself has a finite 
size. Therefore, the anisotropy due to bubble growth might not be 
very important. 

In Fig. 3, we see that the quasar H II region in the column 
density map looks more distinct from other H II regions than in 
the thin (0.64 cMpc for the 164 cMpc box and 1.2 cMpc for the 
607 cMpc box) ionization fraction slices in Fig. 2 in the sense that 
the difference in size is more visible. This points to a more equal 
extend of the H II region in all three spatial directions. 

To test the 3 -dimensional structure further, we show smoothed 
differential brightness temperature maps, see Fig. 4. Here we con- 
centrate on the large box case since this comprises the most lumi- 
nous of our quasars and it does not by construction exclude other 
large scale H II regions, see the comment in Sect. 2. To produce 
those maps we assume that the IGM is heated well above the CMB 
temperature so that we can use Eq. 4 to convert the neutral gas den- 
sity into a differential brightness temperature 



5Tb — Cniii{z, comoving) Vz + 1 



(4) 



with C — 4.6 X 10^ K cm^. We smoothed the field with a gaussian 
of beam-size roughly 3 arcmin (2.75 for z = 8.397 and 2.5 for 
z = 7.57) and subtracted the mean to mimic the effect of a missing 
absolute calibration typical for an interferometer. The width of the 
beam corresponds roughly to the LOFAR core resolution at the cor- 
responding wavelength A '-^ 1.97 m and A ^-^ 1.8 m, respectively. 
If LOFAR had a larger collecting area and thus a higher sensitivity, 
it would be possible to get maps like this. 

In Fig. 4, we show three different planes through the quasar. 
To assess the structure in three dimensions better, we concentrate 
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Figure 4. Brightness temperature plots for the large simulation volume. 
Shown are three slices (xy, zx, and yz) 23.0 Mpc after the quasar turn on 
and one slice without the quasar included (upper left corner). The quasar 
H II region is marked with "A". A large non-quasar region is marked with 
"B", for details see text. 



on two sub-regions, the quasar region A and the large H II region 
marked as B. Region B is not a quasar H II region but formed in 
a highly biased region with many massive halos. Since there is no 
central source that dominates the H II region, we first find the cen- 
tre cell of this region by minimizing the sum of brightness temper- 
atures in a sphere with the same radius as the quasar H II region. 

In Fig. 5, we show slices perpendicular to the three principal 
axes for both regions indicated in Fig. 4. Panel (a) shows the quasar 
region A and panel (b) shows region B. The circles have the radii 
found for the respective regions with the matched filter method, 
see Sect. 6. It is visible that the radius found by the matched filter 
method allows for some emission inside the bubble. The size found 
for region B is larger than for region A, Rh{B) ^ 28.7 cMpc 
compared to Rh{A) ^ 24.9 cMpc. However, as will be discussed 
in Sect. 6, the peak in the signal-to-noise ratio is not as distinct, 
see panel (c) in Fig. 10, which is perhaps obvious given the less 
spherical shape of region B as visible in Fig. 5. To quantify the 
degree of sphericity, we use a method introduced by Iliev et al. 
(2008a). Here, rays are cast from the central source and for each 
ray, the distance where the optical depth r along this ray surpasses a 
certain limit Ts is recorded. In Fig. 6, we show histograms for these 
distances for region A with and without a quasar, as well as for 
region B. As the central point for region B, we chose the same point 
as for the matched filter analysis. For each point, we cast 10 000 
rays and use the threshold optical depth Ts = 4.6. We also tested the 
smaller threshold Ts — 1. We found the curves to be very similar, 
confirming the visual impression that the edges of the H II regions 
are sharp. Comparing the H II region around the most massive halo 
(A) with and without quasar, we see that that the quasar H II regions 
is a bit more spherical and larger in size. The non quasar H II region 
B is larger than the quasar H II region but much less spherical. 



Figure 5. Brightness temperature maps as in Fig. 4, colour coded as indi- 
cated by the colour bar. Shown are three slices perpendicular to the princi- 
ple axes for each H II region. Panel (a) shows the quasar H II region, region 
A, panel (b) shows region B. The side length of each slice is 43.9 cMpc. 
The circles indicate the size of the H II region found by the matched filter 
method, see Sect. 6. Note that region B is not in the same planes as shown 
in Fig. 4 
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Figure 6. Histogram of H II region sizes found by casting rays from the 
centre and registering where the optical depth along the ray surpasses Tg. 
Shown are results for region A with (red) and without (green) the quasar, 
and region B (blue). The threshold optical depth for all cases is Ts = 4.6. 
The dashed lines indicate the medians. 



3.3 Helium ionization region around the quasar 

In Fig. 7 we show the He II fraction for the late quasar turn on. We 
chose this case and not the large box, since it has a higher resolu- 
tion. The circle indicates the calculated He III radius of the quasar, 
R ^ 19.1 cMpc, assuming all photons able to double ionize helium 



(^7'>54.4ev = ^'^ X 10 ) ^o SO. Sincc the fraction of photons 
going into the ionization of a species depends on its contribution to 
the optical depth, this is an overestimation. 

There are mainly two differences between the H II region and 
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Figure 7. He II fraction for the late quasar turn on case, colour coded as 
indicated by the colour scale. The circle indicates the radius calculated on 
the basis of photon counts, see text. 
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Figure 8. Comparing of the signal in the visibility from a single quasar H 
II region and 1200 hours LOFAR system noise. Green circles represent the 
noise in the real part of the visibility as a function of baseline U for 1200 
hours of observations in a single frequency channel of width 38 KHz at 
redshift 7.57. Red crosses are 20 times the expected signal from a single H 
II region of radius 20 cMpc embedded in a uniform H I density field with 
xm = 0.55 at redshift 7.57. 



the He III region: 

(1) Since the contribution to photons with energies above 54.4 eV 
(the threshold energy of helium double ionization) of stellar sources 
is negligible, the quasar dominates the double-helium ionizing pho- 
ton output. There are no other large-scale He III regions to merge 
with. Therefore, the He III region is much more spherical than the 
H II region. 

(2) The He II - He III front is much shallower than the H II region 
front. 

As discussed by McQuinn & Switzer (2009) and Bagla & 
Loeb (2009), the hyperfine transition line of ^He II at 8.7 GHz 
could be in principle used to detect singly ionized helium. How- 
ever, as McQuinn & Switzer (2009) calculate, the spin temperature 
of the transition is not significantly decoupled from the CMB tem- 
perature. Therefore, this line will only be observable in absorption 
against, for example, another background quasar (Bagla & Loeb 
2009). 



4 SIMULATING VISIBILITIES 

The quantity measured in radio-interferometric observations is the 
visibility V{U,iy) which is related to the specific intensity pattern 
on the sky 1^(0) as 



/ 



V{U,iy) = / d^OA{0)h{0)e 



(5) 



Here the baseline U = d/X denotes the antenna separation d pro- 
jected in the plane perpendicular to the line of sight in units of the 
observing wavelength A corresponding to the observing frequency 
jy, 6 is a two dimensional vector in the plane of the sky with origin 
at the centre of the field of view (FoV), and A{0) is the beam pat- 
tern of the individual antenna. For LOFAR, we approximate A{6) 
by a top-hat function with FoV 5°. 

The visibility recorded in radio-interferometric observations 
of a H II region during the EoR is actually a combination of several 
contributions 

V{U, ly) = S{U, ly) + HF{U, u) + F{U, u) + N{U, z/), (6) 



where S(U, f) is the signal from an H II region of comoving radius 
Rh embedded in an uniform H I distribution. HF{U ^ u) is the con- 
tribution from the fluctuating H I outside the H II region which arise 
because of the existence of other H II regions and density fluctua- 
tions in the H I. F(U, v) and N(IJ ^v) are the contributions from 
the foreground and the system noise, respectively. The H I fluctua- 
tions cannot be reduced by increasing the observing time and thus 
put a fundamental lower limit on the size of H II regions that can 
be detected (Datta et al. 2007, 2008). The contribution from the 
fluctuating H I is less than the signal if the targeted H II region is 
bigger than the average size of surrounding H II regions (Datta et al. 
2008). Foregrounds F(U ^v) are much stronger (several orders of 
magnitude) than the signal, but they are expected to be smooth in 
frequency. In Sect. 7 we show that foregrounds can be subtracted 
out and the residual foreground is below the signal for LOFAR (also 
see Datta et al. 2008, for the GMRT and MWA). In the next section 
we discuss the signal in visibility and how we calculate it from our 
simulated image cubes. 



4.1 Simulating Signal 

The expected visibility pattern for a simplistic case such as a spher- 
ical H II region embedded in an uniform H I density can be calcu- 
lated analytically (Datta et al. 2007). A single frequency channel 
will observe a circular ionized region across the H II region. The 
visibility for any circular ionized region of radius Rj, at a comov- 
ing distance r^ can be written as. 



S{U,iy) 



-.Lxn,eW-^-'^ [^4^^^] e (i - ^) (7) 



where I^ is the mean specific intensity for the mass average neutral 
fraction xhi = 1, ^i/ = RvJTtj is the angular size of the circular 
ionized region across the H II region. The phase factor g^TriC/.^c 
arises because of arbitrary position of H II region Q^ with respect 
to the centre of the FoV J\(x) and G(x) are the first order Bessel 
function and Heaviside step function, respectively. The radius R^ 
of any circular region through the H II region is calculated from 
the H II region radius i?b by R^, — Rh ^l — (Az//Az/b)^ where 
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Figure 9. LOFAR coordinates of the 48 core stations (panel a). The LOFAR uv tracks for a 4 hours observations towards the NCP (panel b). 



Az/ = i/ — z/c is the distance of the circular region from the H II 
region centre Vc and Az^b = i?b/^ is the H II region size in 
frequency units. The U range, frequency range and peak value of 
a signal S{U ^ v) scale as Rh~^, Rh and R^, respectively i.e, for 
a larger H II region the signal is confined to shorter baselines. Fig. 
8 plots the expected signal S{U,iy) multiplied by 20 as a function 
of baseline U for the central circular ionized map of a single H 
II region of radius 20 cMpc embedded in a uniform H I density 
with xhi = 0.55 at redshift 7.57. For further details we refer the 
reader to figures 1 & 2 and the discussion in the section 2.1 in 
Datta et al. (2007). In our current simulations described in Sect. 3 
we see that the quasar H II region is not spherical and the outside 
H I is fluctuating. However the discussion above is still useful to 
understand how the signal looks as a function of baseline U and 
frequency u. In addition, the filter we consider later is based on the 
signal we discussed above. 

The radiative transfer simulations provide us with the bright- 
ness temperature distribution 6Th{0, v) at different angular loca- 
tions and frequencies u assuming a global heating as described in 
Sect. 3. The 'early' and 'late' simulation cubes which have a side 
length 163 cMpc are ^^ 1° and ^ 10 MHz wide in angle and fre- 
quency, respectively. For the large box simulation, these numbers 
are ^ 4° and ^ 36 MHz. The image 6Th{0, v) at each frequency 
channel is multiplied by ^^ = ^^^2^ to convert to the specific 
intensity I^ (0). Next we perform the 2D discrete Fourier transform 
(DFT) of the product to simulate the complex visibilities. The LO- 
FAR core has baselines in the range 20 < U < 2000 which is 
adequate to capture the H I signal from H II regions which are ex- 
pected to be confined to small baselines U < 1000 for the stages 
of reionization we consider here. 



4.2 Simulating Noise 

In this section we describe our method of simulating the expected 
noise for LOFAR in the visibility. 

The system noise contribution N{U, v) (see Eq. 6) in each 
baseline and frequency channel is expected to be an independent 
Gaussian random variable with zero mean {{N) = 0) and whose 



rms is independent of U . The predicted rms noise contribution in 
the real (or imaginary) part of a single visibility for a single polar- 
ization can be written as (Thompson et al. 2001) 



■'- ' rms 



V2r]kBTsys 

AeffVAzyAt 



(8) 



where Tsys is the total system temperature, /cb is the Boltzmann 
constant and 77 is the antenna efficiency. We assume 77 to be 1 . Aes 
is the effective collecting area of each antenna, Au is the frequency 
channel width and At is the correlator integration time. Eq. 8 can 
be rewritten as 



7V^ 

1 y ri 



i-95Jy(ii)(^"(^)-^^^(f) 



-1/2 



(9) 



For LOFAR we assume 



Tsys = Tsky + Treceiver = [60 {u /300M H z)'''-^^ + 140] K 

(Harker et al. 2009) and Aes = 526 (150 MHz/z/)^ m^ ^ For two 
polarizations the noise is reduced by a factor of V2. Our small and 
large simulated cubes have resolutions 0.64 cMpc and 1.2cMpc 
which are equivalent to frequency resolutions Av ^ 38 KHz and 
^ 71 KHz respectively at redshift z = 7.57. For Au = 70 KHz, 
integration time At = 100 s, this gives A^r^s = 608 mJy for two 
polarizations at redshift z = 7.57. 

The baselines obtained using DFT of our simulations are uni- 
formly distributed on a uv plane. In real observations, the baselines 
will have a complicated distribution depending on the antenna lay- 
out and the direction of observations. We use 48 high band antennae 
(left panel of Fig. 9) of the LOFAR core to simulate the uv tracks 
for 4 hours of observations with 100s integration time at declina- 
tion 90° (NCP^) (see right panel of Fig. 9). We then count the total 
number of measurements n{U) in the area defined by u±Au/2 and 
V ± A^;/2 around the baseline U{u, v) obtained using DFT of our 

^ http://astron.nl/radio-observatory/astronomers/lofar-imaging- 

capabilities-sensitivity/sensitivity-lofar-array/sensiti 

^ One of the potential fields for the LOFAR -EoR 
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simulations, where Ai^A^; is the grid size determined by the simu- 
lated cube size. We then reduce the noise rms in both the real and 
imaginary part to N^^/ J n{U). We generate a Gaussian random 

number with rms N^^/ J n{U) for each grid point. We assume the 
signal and noise to be zero where there are no measurements. We 
do not expect the rms to change much for different choices of inte- 
gration time. Since we show our results for a total observing time 
of 1200 hours, we further reduce the noise by a factor of V300. 

Finally we simulate 1000 independent noise realizations to in- 
vestigate how the system noise changes the H II region size deter- 
mination and quantify statistically the accuracy at which different 
H II region sizes can be measured. We discuss the effect of system 
noise on the size determination of a quasar H II region in Sect. 6.2. 



5 MATCHED FILTER FORMALISM: A BRIEF 
SUMMARY 

The details of the matched filter technique were discussed in Datta 
et al. (2007, 2008). Here we briefly describe the main features rel- 
evant for our study. The signal component 5'([/, z^) in the observed 
visibilities F ([/, z^) is expected to be buried deep in other contribu- 
tions (see Eq. 6), some of which are orders of magnitude larger. In 
order to enhance the signal to noise ratio (SNR) we sum the entire 
observed visibility signal weighed with the filter. The estimator E 
is defined as. 



Y,^^dUa,iyb)V{Ua 



T^b 



/E 



(10) 



where S{{U,iy) is a complex filter which has been constructed to 
detect a particular H II region and the '*' denotes the complex con- 
jugate, V{Ua,i^b) refers to the observed visibilities and Ua and ut 
refer to the different baselines and frequency channels in the obser- 
vation. The filter 5'f ([/, z/) is a function of the comoving radius i?f, 
redshift Zi and angular position O^. We do not show this explicitly, 
the values of these parameters will be clear from the context. 

Using the simulated visibilities, we evaluate the estimator as 



E = (AuAvAiy) ^ S;{U^, VJ)V{U^, zy,)/>iv([7^, Vj 



(11) 



where the sum is now over the baselines and frequency channels in 
the simulation. Au, Av, Av are the grid sizes along u,v and v axis 
obtained from the simulation. The product AuAvAi/pN{U, z^) is 
the fraction of all baselines in the interval U, U -\- dU and zy, 
V -\- dv. Note that the normalized baseline distribution function 
Pn{U, z/) is usually frequency dependent, and it is normalized so 
that f d^UdiypN{U,T^) = 1. 



The filter 5'f ([/, z/) is defined as 



(12) 



where S{{U, v) is the expected signal of the H II region of radius 
i?f , at redshift Zi and spatial position 0^. According to the matched 
filter formalism, the SNR is maximized when the filter parameters 
match exactly the H II region parameters i.e, Ri = Rh, Zi — z\, 
and di — Oh. Therefore, the peak in the SNR over Rf curve corre- 
sponds to H II region radius R^. We must note that the above filter 
is slightly different from the one we propose in Datta et al. (2007). 
In Datta et al. (2007) (see Filter I), an extra term is introduced to 
subtract out the smooth foreground components. We drop this term 



here but re-include it later when we discuss the foreground subtrac- 
tion in Sect. 7. We also drop the overall factor {v/vcY from the 
filter. This term accounts for the fact that pn{U ^v) changes with 
frequency. As will be demonstrated in Sect. 7, the residuals from 
the foreground subtraction are below the signal for LOFAR even if 
we drop this term from the filter. 

The variance {{AE^) of the estimator is a sum of three con- 
tributions 



((A^)^) = ((A^)^)hf + ((A^)')fg + ((A^)'). 



(13) 



where ((A^)^)hf, ((A^)^)fg and ((A^)^)n are the contribu- 
tions from H I fluctuations, foreground residual and the system 
noise, respectively. To calculate the variance from the H I fluctu- 
ations ((A^)^)hf one would need to have many different real- 
izations of the neutral density field surrounding always the same 
H II region. Simulating many EoR simulations of the kind we 
present here is computationally expensive and hence we do not 
discuss this here (see Datta et al. (2008) for approximate values 
of ((A^)^)hf). Since the quasar H II region is much bigger than 
the average size of surrounding H II regions, the contribution from 
((A^)^)hf) is expected to be much smaller than the signal {E) 
(Datta et al. 2008). The system noise variance {{AEY)n can be 
calculated analytically using the relation 



{(^Ef)i 



''I''"! 



diypN{U,iy)\Si{U,iy)Y 



(14) 



where (^'^ — \ \ {^^ ) /^a,b. To test our noise simulations de- 
scribed in Sect. 4.2, we calculate ((A^)^)n from 1000 indepen- 
dent noise realizations and also using Eq. 14 for the same set of 
parameters and we find that they agree with each other within 5%. 
This test gives us the confidence that our noise simulations are cor- 
rect. We use the simulated noise to further predict the constraints 
LOFAR should be able to put on different H II region sizes. We will 
discuss this in Sect. 6.2. 



6 RESULTS OF THE MATCHED FILTER TECHNIQUE 

Here we discuss the detectability of the quasar H II regions in our 
simulations for LOFAR. Next we discuss if LOFAR can identify 
quasar H II regions blindly. We discuss the effect of system noise 
on the size determination of an H II region and the accuracy in size 
determination we can achieve for LOFAR. 



6.1 Detectability 



Fig. 10 shows SNR= {E)/yJ {{AE)'^)^ as a function of filter size 
Ri for our different simulations for 1200 hours of LOFAR obser- 
vations. Note that here we do not add system noise to the signal 

but calculate the expected noise rms ^/((A^)^)^ using Eq. 14. In 
Sect. 6.2 we discuss in detail the effect of noise when we add noise 
to the signal. Here we take the filter position (^f , z^f ) to be the same 
as the H II region position (^b, ^b)- This corresponds to a targeted 
observation i.e, the central quasar's position and redshift are known 
from an optical/near infra-red survey and we point the radio tele- 
scope to it (see Datta et al. (2008) for blind search). The bottom 
line (red) in the right panel (c) shows the result for the quasar H II 
region in the large box simulation. We see that the SNR peaks at 
filter size Ri = 24.6 cMpc and the peak SNR is 4.6. According to 
the matched filter principle the SNR peaks when the filter matches 
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Figure 10. SNR={^}/a/{(A^)2)^ as a function of filter size R^ for 1200 hours of LOFAR observations. Panel (a): the two curves (from top to bottom) 
represent the late simulations without the quasar ( see upper panels of Fig. 2b for the images) at two time steps corresponding to redshifts 7.57 and 7.664 
respectively. Panel (b): same as panel (a) for the late quasar turn-on case with quasar (see lower panels of Fig. 2 b for the images). Panel (c) shows the same 
for the large box simulation. The lower curve shows for the quasar H II region (marked with 'A' in Fig. 4) surrounding the most massive DM halo. The upper 
curves shows the SNR for the adjacent H II region created by the stars only (marked with 'B' in Fig. 4). The two thin (black) lines are for spherical H II 
regions of radii 24.9 cMpc (bottom) and 28 cMpc (top) embedded in a uniform H I density at redshift 7.57. The lines are scaled to compare with the large 
box simulation results. 



^ 2 

C/5 





/ ...- ■X-.. 


■ 


fa) 


;;;■'■' 


'*'v 


. / 


/ ^ 


">K... 


■■■:: ' : 


'i^ 






^^: 








■ 



10 



15 



20 
Rf/Mpc 



25 



30 



35 




25 30 35 
Rf/Mpc 



Figure 11. The effect of adding the system noise to the signal. The left panel (a) shows results for the late simulation with quasar for 1200 hours of LOFAR 
observations. Thick solid line shows SNR for different filter sizes Ri for the signal only. Other lines represent randomly selected cases when we add different 
noise realizations to the signal. Right panel (b) shows the same for the large box simulation with quasar. 



best with the signal. Therefore the above Ri value is interpreted as 
the H II region size found by our method. We also calculate the ex- 
pected H II region size in the simulation based on the total number 
of ionizing photons (see Sect. 3 ) which gives 25.1 cMpc. We see 
that both agree well with each other. The middle panel of Fig. 10 
(b) shows results for the late quasar turn-on case. The two differ- 
ent lines correspond to two different times after quasar turn on. 
The left panel (a) is the same as the middle but without the quasar. 
The peak SNR, H II region radii obtained both from matched filter 
method and from total photon counts for different simulations are 
summarized in Table 3. We see that LOFAR has good prospects 
to detect our quasar H II regions at redshift 7.57 when the global 
ionization fraction is ^ 50%. As the LOFAR FoV is bigger than 
our simulation box size, the chance of having at least one such de- 
tectable quasar H II region is very good. We also see that the H II 
region sizes obtained from the matched filter method are in good 
agreement with the expected H II region sizes from total ionizing 
photon counts. Therefore we conclude that the H II region size ob- 
tained from the matched filter method can be used to calculate the 
total number of ionizing photons within a spherical region of size 
equal to the matched H II region size. This can be further used to 
constrain the quasar ionizing photon rate and age. We discuss this 
in Sect. 8 



Table 3. Summary of results obtained from the matched filter method. 





quasar age 


z 


Peak SNR 


H II region size 
(from filter) 


H II region size 
(from total photon) 


early quasar 


23.0 Myr 


8.397 


1.2 


11.6cMpc 


12.0 cMpc 


late quasar 


11.5 Myr 


7.664 


2.2 


16.0 cMpc 


16.4 cMpc 




23.0 Myr 


7.570 


2.8 


19.4 cMpc 


19.7 cMpc 


large box 


23.0 Myr 


7.570 


4.6 


24.9 cMpc 


25.1 cMpc 



The upper (green) thick curve in the panel (c) in Fig. 10 rep- 
resents the large H II region very close to the quasar H II region 
in our large box simulation (marked with 'B' in Fig. 4, see also 
Fig. 5, panel b). We discussed in Sect. 3 that the H II region is 
bigger than the one around the luminous quasar. We find the H II 
region size to be ^ 28.7 cMpc and the peak SNR achievable by 
LOFAR for 1200 hours of observations is 5.8. Therefore, we argue 
that LOFAR cannot alone distinguish between H II regions made 
by clustered galaxy sources and H II regions made by a luminous 
quasar by measuring the size of the H II region. Another possible 
way could be to use the shape of the SNR vs Ri curve to extract 
the information about the sphericity of an H II region. The SNR vs 
Ri curve will be more peaked for a spherical H II region than an 
a-spherical one. The two thin (black) lines in the panel (c) of Fig. 
10 show the results for spherical H II regions of radii 24.9 cMpc 
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Figure 12. The estimator {E) (solid thick red line) for the quasar late simu- 
lation is compared to the system noise rms ^/{{KEy^ (green dashed 
line) for 1200 hours LOFAR observations and the foreground residual 
^/ {{/\EY)yg (magenta and blue dashed lines) for the bandwidth 6 times 
the H II region size in frequency. The magenta and blue dashed lines show 
for the cases when p(t/, u) is frequency dependent and independent respec- 
tively (see text for details). 



(bottom thin line) and 28 cMpc (top thin line) embedded in a uni- 
form H I density at redshift 7.57. The lines are scaled to compare 
with the large box simulation H II regions around quasar and stellar 
sources (marked with A and B in Fig. 4, respectively). As expected 
the SNR vs Ri curves for spherical H II regions are more peaked 
than the simulation. But we note that the quasar H II region is more 
spherical than the stellar H II region. So in principle the shape of 
the SNR vs R^ carries information about the shape of the H II re- 
gions. However, as described below in Sect. 6.2, the system noise 
will make it impossible to use the shape of the SNR-i?f curve to 
make the statement about the sphericity of the regions found. 

Next we discuss the effect of the size of the smaller box that 
we use here. The size of the early and late simulation box is ^ 1° at 
the redshifts of our interest. This is much smaller than the LOFAR 
FoV (^ 5°). The minimum baseline we sample for the early and 
late simulations is t/min ^ 56 for this smaller box. A substantial 
amount of signal comes from the large scales for the H II regions 
we are interested in. In addition, LOFAR has many small baselines 
(see right panel of Fig. 9), so we lose a considerable amount of 
the signal in the estimator. The amount of loss depends on the H 
II region size: for larger H II regions the loss is higher than for 
smaller. We use our analytical equations to calculate the loss in the 
peak SNR for different H II region sizes. For H II regions of radius 
20, 30 and 40 cMpc we lose - 3%, - 4% and - 7% in the SNR, 
respectively. For the large box this will be even smaller and hence 
we do not discuss it here. 



6.2 Effect of System noise on the size determination 

In the previous section we discussed results without explicitly 
adding system noise to the simulated signal, but only calcu- 
lated the expected noise rms in the estimator. Here, we add 
the simulated noise realizations to the signal and then perform 
the matched filtering. The left panel (a) of Fig. 11 shows the 

SNR= {E + Ai;N)/\/((A^)2)N vs Ri for four randomly 
selected noise realizations for 1200 hours of LOFAR observations 



for the late quasar simulation. The solid thick (red) line shows 
the result without the noise. We see that in some cases there are 
clear peaks but shifted to either side of the original peak. We also 
see various features arising because of the noise such as a nearly 
flat SNR curve (no obvious peak), double peaked, suppressed and 
enhanced SNR. In real observations one of these SNR curves 
will be observed. Now we only consider the realizations where 
the peak SNR is higher than 3 (because a SNR less than 3 is not 
considered a detection at all), and then find the best matched filter 
size Rf from the peak SNR. We get an average H II region size 
{Ri) = 19cMpc and gr^ = 4.06 cMpc. The right panel shows 
the results for the large box simulation. In the simulation the H 
II region is detectable with SNR=4.6. As expected, the peaks are 
more prominent in almost all cases and therefore the H II region 
size can be detected with higher accuracy with (i?f ) = 25.3 cMpc 
and (tr^ — 4.0cMpc. H II region size can be measured more 
accurately for higher SNR measurements. We also find that 
(Tr^ — 2.5, 1.75, 1.2 cMpc for our late quasar simulation where 
the H II region radius is 19cMpc if we have 5cr, 7cr and lOcr 
detection respectively. 



7 FOREGROUND SUBTRACTION AND ITS EFFECT ON 
THE DETECTION OF A H II REGION 

Astrophysical foregrounds are expected to be several orders of 
magnitude stronger than the H I signal (Ali et al. 2008; Bernardi 
et al. 2009; Paciga et al. 2011; Santos et al. 2005; Jelic et al. 2008) 
at frequencies relevant for EoR studies. Major components are the 
diffuse synchrotron and free free emission from our own Galaxy, 
extra galactic point sources etc. (see Jelic et al. 2008, for more 
detail). Various methods have been developed to subtract the fore- 
grounds (Wang et al. 2006; Jelic et al. 2008; Harker et al. 2009; 
Bowman et al. 2009; Liu et al. 2009; Ghosh et al. 2011a,b; Mao 
2012; Chapman et al. 2012)). All the methods rely on the fact that 
the foreground components are smooth in frequency. As we first 
proposed in Datta et al. (2007), we modify the matched filter to 
subtract the foreground below the signal. The modified filter is. 



1 r 

BW / 



S{U,v)dv (15) 

i^f-BW/2 

where BW is the bandwidth over which the foreground subtrac- 
tion is done. Vf is the frequency centre of the filter. Note that 
/.7-lw/2^'^f^(^'^)^^ = and Sf'itj.iy) is also symmetric 
around the frequency centre i/f i.e, Sf^ ([/, z/ — z/f ) = Sf^ {U,—{iy— 
z/f)). This means, the filter will subtract any frequency independent 
component in the frequency range z/f — BW/2 to z/f + BW/2. 
Because of the symmetry, the filter will also subtract out any com- 
ponent which is an odd function of (z/ — z/f ) i.e, components pro- 
portional to (ly — z/f) (linear), (u — z/f)^ and so on. However the 
quadratic terms will remain as the residuals. 

We next calculate the contribution of foreground residuals to 
the variance of the estimator using the the analytical formula (see 
Datta et al. (2007) for more details) 



{iAEf)FG = fd^ufdiyif 



dv2 



~dT 



dBjy^ 

dT 



xp{U, vi)p{U, zy2)^^f*(t/, iyi)S^{U, V2) 
XC27TU (1^1,1^2) (16) 



© 2012 RAS, MNRAS 001, 1-17 



14 Datta et al. 



where (^^) converts the temperature to the specific intensity and 
C2'kU (^1 , ^2) is the multi-frequency angular power spectrum of all 
foreground components. 

For each component of the foreground the multi-frequency an- 
gular power spectrum is modelled as 



C2'ku{i^1,'1^2) = A 



/loooV 

\2t^U ) 



(17) 



We assume the foreground components to be correlated at two 
different frequencies h>i and z/2 and ignore possible small decor- 
relations. We consider two major foreground components: dif- 
fuse galactic synchrotron emission from our own galaxy and point 
source contribution (clustered and Poisson part). We have adopted 
the parameter values for A, a and ^ from Santos et al. (2005). Un- 
der the assumption of the above foreground modelling Eq. 16 can 
be written as 



((AE)»),o 
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i/f+BW/2 



i/f-BW/2 



di^Mu,i^)\sf{u,i^)\ 



(18) 



where the index 'i' indicates different foreground components and 
fi{U,iy) = (^) p{U,iy) {"^y". Note that here we assumed 
the baseline distribution p{U,iy) to be circularly symmetric. This 
assumption together with the foreground model we assumed in 
Eq. 17 reduces Eq. 16 from four integrals to two and hence the 
evaluation of Eq. 18 is very fast with very good accuracy. The 
terms (^^), p{U, ly), (77)°"* are expected to be smooth in fre- 
quency. Therefore they can be written as (in the frequency range 
ly - BW/2 < zy < zy + BW/2) 



f{U,iy) 



^2r(t/,^f) 



+ [ly - z^f) 



2! 
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(19) 



where the prime denotes derivatives with respect to v. As we 
have discussed before in this section, the modified matched fil- 
ter Sf^{U, ly) removes the constant part /([/, z/f), the linear term 
{ly — h'i)f' and all terms with odd powers. Only the quadratic term 
[ly — ViY 2\^ and other terms with even powers remain as a 
residual. This is the reason why the filter is so successful in remov- 
ing the foreground below the signal. 

In Fig. 12 we plot ^J {{/\Ey)YG for different filter sizes R^ 
(magenta dashed line) for the bandwidth 6 times the filter size Ri 
in frequency. For comparison we also plot the estimator {E) for our 
smaller box simulation with quasar at 2; = 7.57 and the noise rms 

y ((A^)2)iv. We see that y ((A^)2)fg (magenta dashed line) is 
smaller than the estimator for the filter size Ri < 30 cMpc. We fur- 
ther find that a significant amount of the residual foreground comes 
from the fact that the baseline distribution p{U, ly) is frequency 
dependent i.e, the number of baselines in a uv cell changes with 
frequency. Therefore the effective foreground which basically de- 
pends on the product p{U,iy)iy~^ (see Eq. 18) becomes less smooth 
than the original foreground. Since the baseline distribution at dif- 
ferent frequencies is completely known, the change in the number 
of baselines with frequency for a uv cell can be compensated by 
introducing this in the filter. This way the effect of />([/, v) chang- 
ing with frequency can be avoided. The blue dashed line in Fig. 12 



shows y {{AE)'^)fg for the baseline distribution independent of 
frequency. We see that the residual is much below the estimator for 
the case where />([/, v) is frequency independent. 

We note that the modified filter subtracts out some H I signal 
as well. The amount of H I signal it subtracts depends on the size of 
the bandwidth BW over which the foreground subtraction is done. 
The larger the bandwidth, the smaller is the loss of signal. For H II 
regions of radii between 20 cMpc to 25 cMpc at redshift 7.57 the 
loss in SNR is around 24%, 15% and 10% for the BW 4, 6 or 8 
times the filter size Ri in frequency. On the other hand, the level of 
foreground residuals increases for larger bandwidths. The choice 
of optimal bandwidth depends on the size of H II region and the 
foreground model. In our case we find that a bandwidth 6 times the 
radius of H II region in frequency is reasonable. We notice that the 
modified matched filter does not peak at the H II region radius but 
at a slightly large size. We find that the SNR peaks at filter radius 
Ri = 25.4 Mpc for a spherical H II of radius Rb = 25 Mpc at 
redshift 7.57 for a bandwidth 6 times the H II region radius. This is 
a small effect compared to the uncertainty from the system noise. 



8 DISCUSSION AND CONCLUSIONS 

We investigated three cases of a luminous quasar shining during 
the epoch of reionization in the most massive halos of our simu- 
lation volumes. Two cases used a smaller simulation volume with 
a side length of 163 comoving Mpc and one used a substantially 
bigger volume with a side length of 607 comoving Mpc. One of 
the cases in the smaller volume considered a quasar turn on at an 
early stage of reionization, when the global ionization fraction was 
about 25%. The other two cases considered a quasar turn on at a 
later stage of reionization, when the global ionization fraction was 
about 45%. We found that the 163 Mpc simulation volume was es- 
sentially too small to be able to host many large scale H II regions 
at the same time and therefore, we concentrated on the larger sim- 
ulation volume to investigate if the quasar changes the size and the 
morphology of the H II region in a way which makes the quasar 
H II region distinctly different from H II regions solely formed by 
clustered galaxy sources. We found that although the quasar domi- 
nates the ionizing photon output in, and the growth of the H II re- 
gion hosting the quasar, other highly overdense regions with many 
clustered massive halos can produce H II regions similar or even 
larger in size than the H II region hosting the quasar. Although de- 
tailed investigation of the H II regions shows that the extend in all 
3 spatial dimensions of the quasar H II region is more equal (i.e. 
the quasar H II region is more spherical) than large H II regions by 
galaxies alone, given observational constrains, it might not be pos- 
sible to tell the difference using 21cm observations. Furthermore, 
we found that using a quasar spectrum L oc z/~^'^, the sharpness 
of quasar H II regions will not be measurably different from galaxy 
H II regions. 

Next we use the matched filter technique, which combines the 
signal optimally, to investigate the possibility of detecting a single 
QSO HII region in redshifted 21 cm observations using LOFAR. 
We find that with 1200 hrs of observations, LOFAR should be able 
to detect a single QSO HII region of comoving radius 20 Mpc and 
25 Mpc or bigger with at least 3 and 5 a significance, respectively 
at redshift 7.57 when the universe is 50% ionized. We also find 
that the H II region formed by a quasar at earlier stages of reion- 
ization, in our early quasar turn on case, is not large enough to be 
detectable. We calculate the accuracy at which the bubble size can 
be measured in case of a significant detection. For our best case. 
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the quasar in the 607 Mpc simulation, we find for 100 realisations 
of including system noise Rh = (25.3 ± 4.0) Mpc assuming 1200 
hours of LOFAR integration time. 

We also considered the effect of foregrounds. For this, we 
modified the matched filter to subtract the foreground components. 
For a reasonable foreground model the modified matched filter is 
able to subtract the foregrounds such that the residuals are far below 
the signal level. 

Since we found that clustered galaxy sources in highly over- 
dense regions are able to form H II regions of comparable size to 
the quasar H II region, we also used the matched filter technique on 
a non-quasar H II region to investigate if the form of the SNR-curve 
can be used to tell quasar H II regions apart from H II regions solely 
formed by clustered galaxy sources. Although the resulting SNR- 
curves are substantially less peaked, the system noise hides this 
information. We therefore conclude that a blind search for quasars 
using redshifted 21cm maps is not likely to be successful. 

The radii R^ found by the matched filter method for the three 
cases after 23.0 Myr are 11.6 cMpc, 19.4 cMpc and 24.9 cMpc. 
This is in rough agreement with estimates for the expected bubble 
size on the basis of photon counts, R = 12.0 cMpc, 19.7 cMpc and 
25.1 cMpc taking into account both the stellar photons emitted in 
the quasar region and the quasar photons but neglecting recombina- 
tions. We can now use the radii found by the matched filter method 
and Eq.3 to calculate backwards to estimate the rate of ionizing 
photons emitted by the quasar. Since the best prospects of detect- 
ing an H II region by the matched filter method are for the quasar 
in the large simulation volume, we restrict the following analysis to 
this case. 

Without other information on the quasar, we can assume that 
the minimum time ton, the quasar has been shining at the time of 
observation corresponds to the light-travel time from the centre to 
the edge of the H II region, so ton ^ tJJJi^^ = Rh / {1 + z) / c, where 
c is the speed of light. 

Using this lower limit results in an upper limit for the ionizing 
photon rate from the quasar. Rewriting Eq.3 to include an unknown 
contribution of stellar photons A^^ (in terms of total emitted pho- 
tons contributing to the H II region), the ionizing photon rate from 
the quasar can be estimated to 



^7"< [Rli:^nc 



n; 
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c{l + z) 



(20) 



Here, ric is the comoving average number density of atoms in the 
universe. Assuming ric as average number density in the H II region 
is a good estimate for radii larger than about 20 cMpc. 

If we now concentrate on our best case, the 607 cMpc simula- 
tion, we get the following result: 



iV7° < (1.26 X 10^'^ - 3.34 X 10" 



'n;)s- 



(21) 



If we neglect the contribution from the clustered galaxies com- 
pletely, this gives iV^'° < 1.26 x 10^'^ s"\ This is roughly a 
factor 5 more than the ionizing photon rate that was used in the 
simulation. Roughly, a factor 2 comes from neglecting the contri- 
bution from the clustered galaxies, a factor 2.4 comes from the time 
estimate. The la uncertainty from the size estimate translates into 
cr(iV^^°) = 4.0x 10^S-\ 

In case the rate of ionizing photons is known by estimates of 
the total luminosity and assumptions about the unabsorbed con- 
tinuum spectrum (such an estimate was made for the z = 7.085 
quasar by Mortlock et al. 2011, and resulted for that case in an 



ionizing photon rate of Fu 



1.3 X 10^^ s-^), the radii found 



quasar lifetime. Assuming we would be able to find the exact val- 
ues for the photon rates, see Table 2, we would estimate the quasar 
on time to be: 



to 
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(22) 



by the matched filter method can be used to put lower limits on the 



Ignoring the contribution from the galaxy clusters N* results in 49 
Myr. This is roughly a factor of two longer than the actual quasar 
on time, as would have been expected given the total photon num- 
bers from the quasar and the clustered galaxy sources, as presented 
in Table 2. The la uncertainty from the size estimate translates 
into a (ton) = 23 Myr. However, given that current uncertainties 
of quasar ages are higher than this^°, it would still be an improve- 
ment for limiting the quasar lifetime. As Majumdar et al. (2011) 
already point out in their conclusions, the main advantage of using 
the matched filter method is, that one does not need to rely on a few 
line-of-sight measurements but instead, takes the full 3D informa- 
tion into account. 

Obviously our study is limited to the cases we considered. We 
ignore any relativistic effects here, both in simulations and in the 
estimates. The relativistic effects result in smaller H II regions than 
the case where they are not included. It can be quite significant for 
the initial phases of the expansion of an H II region around a lumi- 
nous QSO, but becomes less important as time progresses (Shapiro 
et al. 2006). However, we expect the effect to be small for our cases 
for two reasons: First, the ionizing photon luminosity N^^^ even for 
the brightest QSO in the large box simulation is a factor ^ 4 lower 
than considered (10^^/s) in Shapiro et al. (2006) and secondly, the 
QSO turns on in a big H II region already created by stellar sources. 
Our study is limited to the quasar luminosities we assumed. If the 
quasar luminosity would be substantially higher, the H II region 
it is located in would be substantially larger. This would increase 
their detectability with the matched filter technique and would in- 
crease the relative accuracy of the estimate of their total ionizing 
photon production and life times. More luminous quasars than we 
considered could occur for two main reasons, namely sample vari- 
ance and a higher conversion factor between halo mass and quasar 
luminosity. 

Even though the size of our largest simulation volume ap- 
proaches that of a single LOFAR field, the entire sky presents a 
volume ^ 1000 larger and more massive halos may exist. How- 
ever, such massive halos will live in extremely biased regions, cor- 
responding to high a peaks in the density, with associated high lev- 
els of ionized flux from galaxies. Our results suggest that, for our 
assumed efficiency parameters for both galaxies and quasars, the 
quasar contribution to the total number of photons ionizing that 
region will be around 50%, with a decreasing trend for higher a- 
peaks. Consequently, although the HII region around a more mas- 
sive halo should be larger and therefore easier to detect in the 21cm 
observations, it will be considerably affected by stellar photons, 
which will have to be taken into account when analysing its size. In 
addition, if the trend we see is correct, in a relative sense the quasar 



10 Haehnelt et al. (1998) find quasar lifetimes to be 10^ - 10^ yr from 
matching stellar and quasar luminosity functions at z ~ 3. Haiman & Hui 
(2001) also find quasar lifetimes around 10^ — 10^ yr from clustering anal- 
ysis of dark matter halos, modelling luminosity function and matching with 
clustering of quasars. Lu & Yu (2011) discuss that recent studies of the 
proximity effect using Lya forest spectra seem to be contradictory since 
some imply lifetimes < 10^ yr (e.g. Kirkman & Tytler 2008) while other 
imply lifetimes > a few 10^ yr (e.g. Worseck & Wisotzki 2006). 



© 2012 RAS, MNRAS 001, 1-17 



16 Dattaetal, 



may contribute less to such a region, and therefore have a lesser 
impact on it. 

A higher luminosity efficiency for our halos would be possi- 
ble if black hole growth is substantially more efficient than we as- 
sumed. However, for a reasonable Eddington efficiency A, we are 
already close to the upper limit of the halo mass - black hole mass 
relation which assumes all baryons in the galaxy to be located in 
the bulge. Also, the most massive black hole found at 2; = 7.5 in 
the simulations of (Di Matteo et al. 2012), who followed the black 
hole growth in detail in a volume comparable to our large box, is 
not that much more massive than what we assume. So, it would 
seem to require a strong discrepancy between the joint growth of 
the halo and its black hole in order to obtain quasars that are an 
order of magnitude more luminous than we have assumed. 

The quasar could also have a larger impact if the ionizing pho- 
ton output from galaxies would be lower than we have assumed. 
However, current constraints on reionization already limit the effi- 
ciency to factors of a few around the values assumed here. Much 
lower efficiencies would not only take us further away from the 
electron scattering optical depth as measured by WMAP but also 
result in a completion of reionization after z ^ 6. In fact, since 
we did not take into account gas clumping below the resolution of 
our radiative transfer grid nor (at least in the smaller of the two 
volumes concerned) the absorbing effect of Lyman Limit Systems, 
slightly higher ionizing efficiencies for the galaxies could actually 
be needed for reionization to complete before z ^ 6. In any case, 
one may hope that a detection of the redshifted 21cm signal over a 
substantial redshift range will give us an approximate reionization 
history which could then be used to constrain the ionizing photon 
output from galaxies. 

As we have shown, the largest HII regions at a given epoch do 
not necessarily host a quasar and so the detection of such a region 
in 21cm observations cannot be used as a secure method for finding 
quasars. However, since such a region does point to an overdense 
part of the Universe and they are relatively rare, it would still make 
sense to use such them for follow up deep optical/near-infrared 
searches for quasars since they may indeed contain one. Although 
an exploration of the efficiency of processing redshifted 21cm data 
using a full four-dimensional (position and size) matched filter 
search is beyond the scope of this paper, our results do seem to 
suggest this to be a feasible approach. 
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